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Granular flows down inclined channels with smooth boundaries are common in nature and in- 
dustry. Nevertheless, flat boundaries have been much less investigated than bumpy ones, which 
are used by most experimental and numerical studies to avoid sliding effects. Using DEM numer- 
ical simulations with side walls we recover quantitatively experimental results. At larger angles 
we predict a rich behavior, including granular convection and inverted density profiles suggesting a 
Rayleigh-Benard type of instability. In many aspects flows on a flat base can be seen as flows over an 
effective bumpy base made of the basal rolling layer, giving Bagnold-type profiles in the overburden. 
We have tested a simple viscoplastic rheological model (Nature 2006, vol 441, pp727-730) in average 
form. The transition between the unidirectional and the convective flows is then clearly apparent 
as a discontinuity in the constitutive relation. 

PACS numbers: 47.57.Gc, 45.70.-n 



I. INTRODUCTION 

Granular dense flows down inclined channels preserve 
the complexity of granular flows while remaining simple 
enough for a detailed analysis pQ. They are of interest 
in engineering applications involving conveying of solid 
materials such as minerals, or in geophysical situations 
like rock avalanches or pyroclastic flows. This article fo- 
cuses on flows down a flat, frictional incline. These flows 
differ substantially from those on a rough or bumpy base 
with macroscopic asperities on the order of the diameter 
of the flowing particles. We have developed simulations 
that model the experimental configuration used by Louge 
and Keast j2] : shallow flows (as 7 grain diameters at rest 
height, see Fig. [4]) in a wide (ss 68 grains) chute with 
flat frictional surfaces. We investigate flows in a range 
of inclination angles containing the range of experimen- 
tally observed Steady and Fully Developed (SFD) flows. 
We reproduce the flow properties quantitatively and an- 
alyze the internal flow structure. We show that above 
a given inclination angle granular convection occurs in 
association with inverted density profiles. To our best 
knowledge our work is the first to predict that secondary 
flows also exist with flat boundaries for SFD flows. The 
basal rolling layer can be seen as an effective "bumpy" 
base for the core flow sliding on top of it [H |3] ■ In these 
conditions, we show that velocities in the main bulk of 
the flow follow a Bagnold scaling. This type of flows is 
associated with a constant homogeneous inertial num- 
ber for SFD flows on a bumpy base [I]. This led us to 
study the theology of these flows, and test whether the 
viscoplastic rheology holds. 

The paper is structured as follows. The next section, 
which can be skipped by specialists of the field, is devoted 
to the state of the art. Section 3 gives details about the 
simulation method we will use. Global properties of the 
flows are studied in section 4. Section 5 is devoted to de- 



tailed results concerning the packing fraction, pressure, 
velocity and "granular temperature" fields. The rheolog- 
ical study is presented in Section 6. Concluding remarks 
are given in Section 7. 



II. THE STATE OF THE ART ON GRANULAR 
FLOWS DOWN INCLINED CHANNELS 

Significant progress have been made during the last 
decades in describing dense granular flows, nevertheless 
they continue to resist our understanding and remain an 
active field of research. Very dense quasistatic regimes 
are usually described by plastic models [5], and a kinetic 
theory of granular gas has been developed [6] that can ac- 
curately render the behavior of dilute flows. A viscoplas- 
tic description for dense fluid regimes has been proposed 
based on a dimensionality analysis in the unidirectional 
case [I]. This rheology has then been extended to 3D 
for incompressible flows jj]. In all these cases, when the 
parameters are set according to the expected theoretical 
values {e.g. high packing fraction dense uniform flows, 
initial and boundary conditions, etc.), the proposed con- 
stitutive equations match the experiment {e.g. collapse of 
velocity profiles |8_ ) . Extensions are then proposed to ac- 
count for variations of nearly related cases: an extended 
granular gas theory taking into account correlations in 
denser cases [3], or a variant of the viscoplastic model 
for compressible flows [TU]. Despite all these efforts, a 
comprehensive theory is still missing in the general and 
most common case where the coexistence of both dense 
and dilute parts are observed within the same flow, and 
which would correctly incorporate the influence of bound- 
ary conditions such as side walls and bottom. 

A large corpus of studies exists on dense granular flows 
down an inclined plane chute (more than 100 references 
in [4], additional ones in (T|). The boundaries are known 
to change the flow structure [TT|. The choice of wide 
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channels is an attempt to avoid the influence of sidewalls. 
In the same way, numerical simulations in periodic cells 
attempt to study flows down infinitely long and wide 
chutes. Most experimental and numerical works avoid 
the inherent discontinuity and sliding at the base by cov- 
ering the surface with glued, fixed grains of the same 
nature as these involved in the flow [5J [TTHT3] . In these 
conditions there exist limits on the lower inclination an- 
gle and the piling height below which the grains do not 
flow. Above these thresholds and for moderate inclina- 
tions, dense fluid flows present a negligible velocity at the 
bumpy base. Thin SFD flows comprising a few layers of 
grains exhibit a nearly linear, sheared vertical profile of 
the velocity (g], 2D experiments [T4HT6] . 2D [17] and 3D 
periodic |12| numerical simulations). For thicker flows, a 
Bagnold scaling is observed in the core of the SFD flow, 
with lower velocities at the base ([3], 3D periodic numer- 
ical simulations [12]). At larger angles of inclination the 
flows are more dilute and an inverted density profile is ob- 
served [TUl [TBI [T5] . This inversion was analyzed by means 
of the granular gas theory to induce a Rayleigh-Benard 
type of instability [T31 IT9] . The convection rolls take 
the form of longitudinal stripe patterns |101 [T5] . These 
can be reproduced numerically using Periodic Boundary 
Conditions (PBC) [10] provided the width the periodic 
cell, W, is large enough compared to the grain diame- 
ter D for the convection rolls to appear. Convection has 
never been observed in numerical works using W = 10 D 
[El [20]. Below about W « 50D, there seems to be not 
enough space for developing convection rolls |TU] . 

Beside these studies of "unconfined" flows in large chan- 
nels, extensive measurements highlighting the influence 
of walls were performed |21H24| . Experimental and nu- 
merical studies both in 2D and 3D configurations [23] 
show that frictional lateral walls alter the flow proper- 
ties. For instance, SFD flows on bumpy bases are ob- 
served up to large inclination angles where accelerated 
ones are usually expected. Moreover, at any given in- 
clination angle, there is a critical flow rate above which 
a static heap forms along the base. The heap is stabi- 
lized by the flow atop it [22] . Flows atop this sidewall- 
stabilized heap (SSH) differ from SFD flows on bumpy 
base as they occur over erodible bases, but still display 
SFD features. The effect of side walls on SFD flows on 
top of a static pile in a channel has been studied by car- 
rying out experiments in setup of different widths, up to 
600 particle diameters [23]. They show that these flows 
are entirely controlled by side wall effects. 

The bumpy bases made of glued grains case is thus rel- 
atively well-studied. However most industrial conditions 
involve flat boundaries, as well as natural flows occurring 
on smooth bed rocks at the scale of the grains. Surpris- 
ingly very few studies [2] [TTJ I25H32] have considered this 
more common case of flat frictional surfaces. Early exper- 
imental works mention increased flow velocity and slid- 
ing conditions at the boundaries compared to the bumpy 
walls case [UJ. Differences with the bumpy case situa- 
tion are manifest in the flow properties. Velocity profiles 
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Figure 1 . Top view of free surface velocity profiles for granular 
flows down flat and frictional base (reprinted from [2]). The 
influence of the side walls is apparent over 2/3 of the flow 
width. 



(transversal and in height) involve a slip condition at the 
boundaries [29 , compared to the null velocity condition 
at the interface with the bumpy walls. Some unexplained 
surging waves are occasionally observed at the surface 
El EI], blurring its exact location by a layer of grains in 
saltation. More recently Louge and Keast [2] conducted 
experiments on a flat base with a well documented set of 
parameters, with a more detailed analysis than previous 
experimental works on the topic |27ff30] . They confirmed 
the aforementioned observations regarding the flow struc- 
ture and velocity profiles. A layer of rolling grains with 
intermittent jumps develops on the flat base, with the 
rest of the flow sliding on top of this basal layer. The 
influence of the distant side-walls is negligible in terms of 
induced friction: "the relative contribution of side walls 
in the force balance [. . . ] never exceeds 7%" [2J. However 
this influence of the side walls is clearly apparent over 2/3 
of the flow width (Fig. [TJ , reprinted from [2]), leading to 
the conclusion that some other mechanism is involved for 
a long-range influence of the boundary conditions. 

They also reported a range of inclination angles 
[#min, 0max] for the observation of steady fully devel- 
oped (SFD) flows, independent of the flow height, that 
presents a much lower bound than in the bumpy sur- 
face case [S]. The upper bound for SFD flows is also 
provided, but as correctly pointed out by 23] the attain- 
ment of SFD flows is restricted by the physical length of 
the chute, hence so is the maximal angle above which an 
"accelerated regime" is observed. Numerical simulations 
can be used in order to complement these experimental 
results, but literature on numerical studies over flat fric- 
tional surfaces is sparse. Early simulations are reported 
in [25] (2D) and [3TJ (3D). Given their limited computa- 
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Velocity, m/s Packing Fraction 

Figure 2. . Velocity (left) and packing fraction (right) profiles 
obtained by Walton for moderate flow heights. Solid lines 
are instantaneous profiles, dotted line is time average profiles. 
Reprinted from |31[ . 

tional power, implying the use of a small periodic cell and 
a low number of grains, and given their use of monodis- 
perse grains, a direct comparison with the experiment is 
difficult. More recent numerical works considering a flat 
base [5] do not provide a detailed analysis of its influ- 
ence. Without sidewalls SFD flows on a flat base can 
only be sustained for inclination angles whose tangent is 
less than the friction coefficient [2] [31]. Thus, with PBC, 
the maximal inclination angle # max is fixed by the solid 
friction on the base: fi — tan# max . Walton [3T] got effec- 
tively SFD flows for inclination angles 9 whose tangent is 
smaller than the friction coefficient /x, else they acceler- 
ate unboundedly. Nevertheless the experimental value of 
0max lead to choosing a large value of [i which is not com- 
patible with the friction coefficients measured in impacts 
[5J|33]. The value of the lower bound 9 m i n is not available 
in |31j . Velocity profiles as a function of distance from 
base at low angle in [31] show a seemingly linear shearing 
in the bulk region (above the rolling layer - see Fig.[2| for 
thin flows, which turn into const ant- velocity crystallized 
plugs at larger thickness. 

At larger angles but in 2D [23] , the profiles of the pack- 
ing fraction, the velocity and its fluctuations, are of the 
same type as these predicted by kinetic theories (type 
II in [20]), with an inversion of the density profile and 
a higher "granular temperature" at the base than at the 
surface. The state of art on granular flows down flat fric- 
tional channels thus remain largely incomplete, with lim- 
ited numerical simulations not able to complement and 
detail the inner details of the flows reported in experi- 
mental works. 



III. SIMULATION METHOD 

We perform 3D numerical simulations of granular flows 
using molecular dynamics (MD) that consist in integrat- 
ing the equations of motion over time. Each grain is 
represented by a sphere whose diameter is drawn from 
a uniform distribution around the mean value D. The 




Figure 3. (Color online). Overlapping spheres contact model. 
Grain i moves with a translation velocity Vi, and similarly 
for j. The force exerted by grain i on grain j during contact 
(characterized by the normal vector n 1- ^ ) is decomposed into 
a tangential component Fl~^ 3 and a normal component F^ 3 . 
Both depend on the overlap <5 according to a contact model 
detailed in the main text. 



grain model consists of a non-deformable sphere |3"5] 
of uniform material with density p. Deformations are 
taken into account by the contact model, which links the 
normal force F„ acting on each grain to the overlap 5 
that occurs between the non-deformed spheres when the 
grains centers are closer than their diameters would al- 
low (Fig. [3). The linear visco-elastic approach [3S] is 
used: = (k n S + 7„w„) n l_>J ', with n <_>3 ' the con- 

tact normal (unit vector from sphere centers i to j), 
V n = ( v i — v j) ■ n 1 ^ 3 the normal component of the rel- 
ative translational grain velocities, k n a model spring 
stiffness and j n a model viscosity (i.e. giving a linear 
velocity-dependent force). A similar model is applied in 
the tangential direction: Fl~^ J — (k t s + ■ftVt) t l ^ J , with 
Vft'^j = (vj — Vj) — v n T^^ 3 the tangential component 
(i.e. with some direction within the tangential plane) 
of the impact velocity, kt and 7t a model spring stiff- 
ness and a model viscosity, and s is a bounded version 
M < | Ft | /kt of the sliding displacement J vtdr in the 
tangential plane since contact time To (37]. Coulomb fric- 
tion |F t | < /x|F„| is enforced on the tangential compo- 
nent, with a model coefficient p. Below that threshold 
the value of Ft is given by the above equations. The 
torque acting on a grain is computed as q = — r (F t x n) 
with r the grain radius. Both force and torque are used 
for integrating the equation of motions Y] F = ma and 
^2 q = Id) with m the mass of a grain, a its accelera- 
tion, / its moment of inertia, and uj its angular veloc- 
ity vector. Numerical integration is performed using the 
Velocity- Verlet scheme. This whole approach is repeated 
for grain-wall interactions with a different set of param- 
eters fcf, ->*•», tT. ^ ■ 

Solid mechanics induces relations between these model 
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Grain diameter 


D = 2.968 (mm) 


Grain mass (glass density p g ) 


m = 3.42 ■ 10" 5 (kg) 


Gravity 


g = 9.81 (m/s 2 ) 


Grain/grain normal restitution 


e 33 = 0.972 


Grain/grain tangential restitution 


ef 9 = 0.25 


Grain/wall normal restitution 


er = 0.8 


Grain/wall tangential restitution 


ef w = 0.35 


Grain(glass)/grain friction 


M 99 = 0.33 


Grain (glass) /wall (aluminum) 


M 9™ = 0.596 


Grain/grain spring stiffness 


k 33 = 2 ■ 10 5 (mg/D) 


Grain/ wall k!^ w = k 33 (glass Young modulus = aluminum) 


Integration time step 


dt = 10~ 4 {yfgD) 



Table I. Simulation parameters matching the experimental 
values from [2]. 



parameters. For a normal collision between two grains 
the damped harmonic oscillator defined by the above in- 
teraction model leads to a contact duration t c during 
which S > (half of the first pseudo-period) . Normal rel- 
ative velocities before and after contact are then related 
by a constant coefficient of normal restitution e n that 
sets 7„. Similarly the tangential spring/dashpot model 
defines a coefficient of restitution et- Equating both du- 
ration times leads to a relation 7k t (n 2 + (In e n ) 2 ) = 
2k n (ir 2 + (In et) 2 ), which corrects the 7k t — 2k n rela- 
tion from [12] when e n ^ ej. Thanks to these relations 
the simulation parameters can be concisely given in Ta- 
ble u 

The correspondence of these parameters to physical 
values is subject to a few simplifications. The most dras- 
tic one is the use of a single model friction coefficient /i 
for all cases of static, kinetic and collisional frictions. We 
had to use the static friction coefficient instead of the 
other ones - as usually done in MD simulations |10[ [T2"] 
- in order to reproduce experimental values. Hypothesis 
for this model/experiment discrepancy given in the liter- 
ature are the presence of long lasting contacts [2] or the 
use of the visco-elastic contact model itself [35]. Even 
then, the static friction coefficient is known to be quite 
sensitive to the surface properties and its determination is 
itself a topic of debate. We used \i" = 0.33 as measured 
in our lab between spent glass beads. Lorenz et al. [34 
had to erode grains by circulating them in their experi- 
mental facilities for two hours before their results became 
reproducible. We used the value they give fi 9W = 0.596 
for the grain/ wall contacts in our simulations, together 
with all their normal and tangential restitution measure- 
ments (see Table |l| . These restitution coefficients could 
be refined for binary collisions using precise velocity- 
dependent measures fitted by more complicated models 
|37| . but this would not necessarily give better global re- 
sults given the multiple- and long-lasting- contacts [37| . 
so we stick to the experimental values given by [2] 154"] . 
Similarly the use of a more complicated non-linear con- 
tact laws (e.g. Hertz) has been proposed but was found 




Figure 4. Sketch of the MD simulation in a configuration 
corresponding to experiments [2] where an inclined plane is 
bounded by fiat side walls and base. 

to be no better than a linear model on a global scale [59] . 
The value of the spring stiffness shall however be related 
to material properties. A link to the Young modulus and 
Poisson ratio is possible for Hertzian contacts [30]- For 
linear models we had to rely on an ad-hoc approxima- 
tion [41] that leads to k n = 3 • 10 5 N/m = 3.35 • 10 6 mg/D 
(comparable to k" = 2.10 5 N/m in [55J). In any case we 
checked that our results are not sensitive to the choice 
of k n provided it is given a sufficiently high value, so we 
then used the more classical value k n = 2.10 5 mg/ D for 
faster simulations |12| . 

IV. FLOW CONFIGURATION 

A. General setup 

The simulation setup is designed to model the experi- 
mental setup of Louge and Keast [2 , with minor adapta- 
tions (Fig. |4| . The calculational space is bounded on the 
base and on the side walls by fixed flat frictional planes, 
and it is free on the top surface (see Fig. [4]). PBC are ap- 
plied in the flow (x) direction as we cannot simulate the 
whole system with current computational facilities (we 
use a period of L ~ 20D, similar to [12 ). Initial condi- 
tions model the dropping of a loose assembly of agitated 
grains at a small altitude. These low energy conditions 
are combined with a mass holdup H — 4 compatible with 
the experimental configurations for all SFD flows in [2 
(the mass holdup H = Jq + °° ^$-dz quantifies the amount 
of matter above a unit surface, with v the volume frac- 
tion). 

Preliminary simulations using PBC along y and a small 
periodic cell size were first performed and were able to 
recover the results of Walton [31] . However, the range 
of angles of inclination for which steady and fully devel- 
oped flows are reached does not match the experimental 
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Figure 5. (Color online). Translation kinetic energy vs dis- 
tance traveled by the grains, showing regimes that appear 
accelerated within the chute length experimental limit, but 
reaching steady states at larger distances. 



results of Louge and Keast [2J. For example, in simu- 
lations we obtain # m i n ~ 6 — 7°, which is much lower 
than the experimental value (15.5°). Neither a modifica- 
tion of the material parameters such as the friction and 
the restitution coefficients, the polydispersity of grains, 
nor the introduction of other sources of dissipation like 
rolling friction [36 gives values of m ; n close the exper- 
imental one. On the contrary, the introduction of side 
walls separated by a gap W with identical material prop- 
erties as the base is able to increase 9 m i n to 9 m - m = 14°, 
reasonably close to the experimental value (15.5°). Note 
also that, the use of the low polydispersity value given in 
[SJ (±0.7% of the average grain diameter) leads to crys- 
tallized blocks in the flow. Due to the periodicity in x 
these blocks tend to persist for a long time. Experimen- 
tally the grains are re-circulated and this corresponds to 
averaging over multiple realizations making the presence 
of blocks inappropriate. One way to get rid of these ar- 
tifacts was to increase the polydispersity, up to 10% for 
all the results presented below. 



B. The flow regimes 

We ran simulations for a range of angles from 13° to 
23° containing the range of experimentally observed SFD 
flows Visual investigation of the simulations shows 
that the main bulk of the flows rests on top of a basal 
layer of grains, for which there is a combination of long- 
lasting rolling contacts with the flat base that are in- 
terupted by short rebounds. Fig. [5] presents the evolu- 
tion of the kinetic energies of the flows over the average 
distance traveled by the grains. Louge and Keast [5] ex- 



perimentally observed a range of angles for steady states 
from 15.5° to 20°, established over distances less than 3m. 
This matches our finding (see the vertical line in Fig. [5]) 
that flows at higher angles would indeed appear accel- 
erated within the experimental limits. The flows stop 
below 9 m i n — 14° which reasonably matches the experi- 
mental value (15.5°) as there are extra factors not taken 
into account in the simulation, like the abundant inter- 
particle dust that was reported experimentally [2J 151] . 
which might then block the flow at low velocities. The 
effect of the side walls on 9 m i n mentioned by Louge and 
Keast [5J cannot be attributed only to additional friction, 
which is negligible given the shallow flow height and the 
large distance between the walls as aforementioned and 
noted in |5] . Some unknown mechanism is thus at work, 
which will be the topic of further studies. 

Simulations with sidewalls lead to SFD flows except, 
maybe, for 9 pa 18° (±1°) where relatively large fluctua- 
tions are visible on the kinetic energy in Fig. [5J These 
fluctuations, which persist over time (we checked it up to 
a distance of pa 13500-D), were also reported experimen- 
tally in [2J, although it is not easy to determine whether 
these match our simulations results: the periodic size 
we use in x is a fraction of the oscillation wavelength 
observed experimentally. Hanes and Walton [33 , who 
use a bumpy base for their experiments, also report a 
phase diagram with an oscillating regime delimited by 
fuzzy boundaries, at the junction of two SFD regimes, 
with the same PBC numerical interpretation difficulty. 
These oscillations take place at the transition between 
two SFD regimes which have very different behaviors. 
For 9 <G [14°, 17°] flows are unidirectional and grains from 
ordered layers. These are visible in Fig. [6] as regions of 
higher packing fraction, as well as in Fig. [7^. The free 
surface of the flow (Fig. |7^) is convex, higher in the cen- 
ter than on the borders. For 9 larger than 19°, secondary 
flows develop (Fig. 17 1, breaking the layer structure (Fig. 
IfjJ, except for the basal layer of rolling grains (Fig. JtJd) . 
The free surface is concave (Fig. ^jp). The flow height is 
enlarged, with a correspondingly lower average density, 
but the flow remains shallow (height pa 10D for a width 
of 68D). The rolls are thus quite flat. Secondary flows 
have only been observed in experiments within a bumpy 
channel, e.g. [11 (using D — 0.5mm beads) and [TU] 
(D = 0.4mm). The stationary state in [TT] was however 
reached at a distance compatible with our results, scaled 
by the difference in D: less than 3m at 9 = 23.6°. 

In order to quantify the effect of the geometry of the 
base, we have carry out the same simulations with a 
bumpy base consisting of fixed grains and otherwise the 
same parameters (including flat frictional lateral walls). 
In these conditions the grains flow only above 22°, with 
an average kinetic energy of pa lmgD. Compared to the 
kinetic energy pa lOOmgD in Fig. [5] at the same angles, 
and given the much lower # m ; n bound in the flat case, 
we immediately see that no direct comparison is possi- 
ble between the flat frictional base and the bumpy one: 
the basal layer of rolling grains significantly reduces the 
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dissipation. Experimentally [TT] also noted a much in- 
creased flow velocity on flat surfaces. For 9 > 22° fixing 
grains on the base prevents the rolls but induce a large 
internal agitation instead. This is consistent with the re- 
sults of Borzsonyi et al. [TU] which report rolls only for 
thicker flows compared to our shallow flow configuration. 
Therefore, the convection regime is accessible for smaller 
systems with a flat frictional bottom than with a bumpy 
one. 



C. Influence of the parameters and of the initial 
conditions 

To study the generality of the above reported results 
we carry out an extensive study of the effect of the model 
parameters and of the initial conditions which is summa- 
rized below. 

The transitions between the regimes and their charac- 
teristics depend slightly on the model parameters (e.g. 
friction coefficients, polydispersity of the grains) but the 
general features of the flows seem robust. For instance, 
additional runs with [i gg = 0.4 instead of 0.33 shifted the 
start and end of the unidirectional flows up by 1°. The 
rolls are robust to polydispersity (tested with D ± 20%). 
They also appear when fx gg = fj, gw (customary setting 
in numerical works [12;), provided both are greater than 
about 0.54. 

The initial conditions we use consist of dropping a loose 
assembly of grains at z = 2D with a low initial velocity 
and some jitter, which we designed to be approximately 
what the grains would have experimentally when leaving 
the open gate in [2], The corresponding initial energy is 
much lower than that reached in steady state (see Fig. 
|5|. We checked the final SFD states are robust to vari- 
ations of the initial conditions provided these induce an 
initial energy smaller than the energy of the final SFD 
state. However we have not studied the use of larger 
energies, as we suspect there may be hysteresis effects 
on the final energy levels. Literature for the bumpy base 
case also reports [TS] that specific regimes exist with high 
initial energy conditions, so presumably this might be a 
possibility for the flat frictional base case as well. 

We varied the mass holdup so as to match the range 
of shallow flow configurations in [2]: low H = 1 in- 
duces an early transition to a dilute phase without sec- 
ondary flows. Convection rolls appear between H = 3 
and H = 4. They persist even for much higher mass 
holdups (tested up to H = 20). 

The existence of a minimal angle # m ; n for SFD flows 
requires the presence of enough layers of grains, as the 
rolling basal layer may accelerate indefinitely with insuf- 
ficient frustrations on the grain rotations (the limit case 
being a single grain rolling on an flat inclined plane). 
The investigation of the parameters and initial conditions 
mentioned above shows that the properties of the flows 
are robust. Therefore, the simulations reported here are 
qualitatively representative of many others obtained with 
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Figure 6. (Color online). Height-averaged packing fraction 
(smoothed over ±0.5D in z) for each angle of inclination. 
The layering is clearly apparent at low angles, and disappears 
in the presence of the secondary flows. 

different conditions. 



V. STEADY STATES 

This section analyzes the main features of the two fam- 
ilies of SFD states : unidirectional and layered flows, 
and flows with granular convection. All the figures pre- 
sented below report time averaged quantities (e.g. ve- 
locity, packing fraction) computed over 500 \f~D~Jg time 
units in steady state, as well as over the periodic cell in 
the flow direction. These results are stable with respect 
to the particular random seed we use between different 
runs. 



A. Packing fraction 

Values of the packing fraction at the base of the flow 
(Fig. [6]) are compatible with the experiments (see for 
example Fig. 6 of [2 , 6 = 16°, 18°, 20° and H = 4). For 
the unidirectional flows Fig. [6] shows that the average 
volume fraction v sa 0.59 does not vary much in z above 
the basal layer, with clear variations around that aver- 
age at each structured layer. The structuration in layers 
of constant packing fraction (above the basal layer) was 
observed by the previous numerical study [3T] with PBC 
along y. We confirm that structure persists with a poly- 
dispersity of 10%, that it was not an artifact of the use 
of single-sized spheres in [3T|. Another difference with 
[3"T] is the presence of 2 lateral layers (Fig. |7| at the wall 
boundaries, inducing some nearby structuration. 

In the convective regime we observe an inverted density 
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Figure 7. (Color online). Packing fraction map, averaged over x. The layered structure at low angles disappears when granular 
convection takes place. 
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Table II. Global effective friction coefficient on the wall (to 
compare with the microscopic jj, gw = 0.596). 



profile (Fig. RjJ, similar to these reported in the literature 
for bumpy boundaries |10( H?| , while the layer in contact 
with the flat boundaries remains clearly separated (Fig. 
§ and [7j)). 



B. Pressure and effective friction 

We did not implement a full calculation of the stress 
tensor, but stresses on the boundaries are easily deduced 
from the forces exerted by the grains during each contact, 
averaged over a 500 \J D/g time window. Let us denote by 
f 9 ^ w (z) the stress vector exerted on one wall, at height z, 
by the grains. The force exerted by one wall on the grain 
flow is simply: F w ^ 3 = - f^ =Q f™ t 9 ^ w (z)dxdz = 
L f^ Q i w ~* 9 (z)dz. Using the subscript n and t to dis- 
tinguish the normal and the tangential components, the 
local effective friction coefficient on the wall is then com- 
puted as fiw( z ) = \\ft( z )\\ I l|fn( z )ll- The global effective 
friction coefficient is computed as ~(Zw = ||Ft|| / ||F„||. 
The profiles in z and the values of the effective friction 
on the walls are shown in Fig. [9] and Table |TT] The ob- 
served friction weakens with depth and is similar to that 
reported in |42j . 

Let us consider a slab at the top of the 
flow (from z and above) as a continuum. 



F n , F n = normal reaction of the wall / base 



w 



Ft , Ft = tangential reaction of the wall / base 




Figure 8. (Color online). External forces balance on the grain 
flow. 



The balance of external forces along z im- 
P lies: ( Jx=o J^o XT v ( x ' y -> z ') dx dydz} pg cos 6 = 
P(z)LW + 2L t™^ 9 {z').e 2 dz' with P(z) the average 
pressure on an xy plane section computed at height z. 
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Figure 9. (Color online). Friction coefficient profile on the 
lateral walls for several value of the inclination angle 9. 



The contribution of the walls to this balance is always 
very small compared to the other terms (of the order of 
one per cent). The vertical profiles of P(z) are shown in 
Fig. [10] for each angle 6, normalized so the value at the 
base equates the mass holdup H we use. They are linear 
over the most part, except for the basal rolling layer and 
the very diluted top consisting of a few grains in ballistic 
motion. Thus, over the main bulk of the flow, the 
approximation P(z) « (H p — z) vp g g cos is excellent 
(with v the average packing fraction on the bulk). The 
corresponding effective flow heights H p are shown in 
Fig. 10 (inset). They confirm the general dilation of 



the flow with the angle 0, matching the general packing 
fraction decrease of Fig. [6] The total frictional influence 
of the walls on the flow can also be quantified. The 
ratio between the weight of the grains and the friction 
force on the walls: fiw H p /W \22\ 143] is at most 0.06 
for 8 = 23°. This justifies the arguments developed in 
[5] and used here for neglecting the friction on the walls, 
for a shallow flow with a large width. 

The pressure at any position on the flat base Pb(u) = 
i^ b (y).e 2 is shown in Fig. 11 normalized so the average 
value is the mass holdup H = 4. In the unidirectional 
regimes grains do not deviate much in y from their tra- 
jectories along x, leading to large pressure fluctuations on 
the base. In the convective regime, the grain circulation 
due to the secondary flows smooth out these differences. 
In each case pressure is maximal at the center of the flow, 
and decreases in the lateral parts near the walls. 



C. Mean velocity and fluctuations 

A continuous mean velocity field v is com- 
puted using the definition by Serero et al. |4"4"] 
for polydisperse systems, and averaged over time: 




i i i r 
0.5 1 1.5 2 

P(z)/(pDgcos6) 

Figure 10. (Color online). Main diagram: Vertical profiles 
of the averaged pressure for the hydrostatic approximation in 
the center part of the curves. Inset: The effective heights of 
the flows, with respect to the hydrostatic approximation H p 
and maximal flow velocity H v . See the main text. 
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Figure 11. (Color online). Normalized pressure computed 
from each contact with the base, smoothed over ±1D in y 



v ( x ) = (E*Li viM** ~ x )) / (Y,j=i k j ( x i ~ x ) 

where x = (ar, y, z) is the 3D position at which to com- 
pute the average velocity, Xj is the position of the cen- 
ter of grain i, and fej is a kernel that distributes the 
mass rrii of grain i over space. We use the uniform den- 
sity kernel fc,(xj — x) = p when ||xj — x|| < r, with n 
the radius of grain i, and elsewhere. We then define 
a "granular temperature" |H] from the velocity fluctua- 
tions: T — | (j|v|| 2 — ||v|| 2 ^ where the overline denotes 
the above weighted averaging. 
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Figure 12. (Color online). Transverse velocity profile, 
smoothed over ±2.5D in Y, averaged over z. 




5 10 15 20 25 

Vx (units of -/gD) in the W/3 central part 

Figure 13. (Color online). Height profile of the velocity along 
x averaged over the central part of the flow. 



For the unidirectional flows Fig. [12] reproduces the 
experimental mean velocity transverse profile (Fig. [TJ, 
where the shearing layer induced by the walls extends to 
about 1/3 within the flow. The shape of the velocity pro- 
file in the unidirectional regime, considering the average 
packing fraction is constant in z (Fig. [61, is compara- 
ble to the experimental measures in Fig. 3 (4) of [25] at 
similar angles. Within the central part we obtain veloc- 
ity profiles (Fig. 13) similar to these given by PBC |31| . 
The velocity reaches a maximal value at height H v then 
decreases rapidly in the sparse ballistic layer of grains. 
The values of H v can be compared to H p in the inset 
of Fig. [10] The 3D velocity profile of the flow can be 
inferred from Figs. [12] and [13] as a faster region in the 
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Figure 14. (Color online). Main plot: Bagnold velocity profile 
collapse. Inset: effective "sliding" velocities of the main bulk 
of the flow with respect to the basal layer of rolling grains. 



center part, sheared vertically, on top of a basal rolling 
layer of grains. The transverse velocity profile (Fig. 12 1 



is sheared through the whole width in the presence of 
secondary rolls. These convey grains up to the center of 



the flow, as seen in Fig. 17 



The basal layer of rolling and bumping grains can be in- 
terpreted as an effective base for the main bulk of the flow 
on top of it. A sliding velocity V s can be defined as the ve- 
locity in the direction of flow at zq = l.hD. corresponding 
to the mean velocity of the grains in the second layer, just 
above the basal grains. Now, let V' x — (V x — V s ) /\fgD. 
z' = (z- z ) /D, and H' = (H p - z ) /D. Bagnold's 
constitutive equation [TJ] states that the pressure P(z) 

relates to the shear rate with oc y/P(z'). Inte- 
grating this relation and assuming that it holds in a 
reference frame moving with velocity V s , we shall have 
VJ. = A(6) x (# /3 / 2 - (H' - z') 3/2 ) , with A{9) a constant 
that is related to the inertial number defined in the next 
section. Fig. [TI] shows that this is indeed the case, that 
all the vertical profiles of the velocity indeed collapse on 
the theoretical curve above the rolling layer. The "sliding 
velocities" V s are shown in inset of Fig. [14] with an ex- 
cellent fit between the measured V s at z = 1.5D and the 
fitted value from the Bagnold profile. V s increase roughly 
linearly with tan 8. 

Globally, flows on flat frictional surfaces can thus be 
decomposed into a rolling basal layer, above which the 
main bulk of the flow follows the classical Bagnold scal- 
ing. This is consistent with the observations reported 
in pQ [3]. Note that, as we have averaged over the trans- 
verse direction, this analysis however tells nothing on the 
internal flow structure visible in the previous sections. 

Fig. [17] shows the color-coded map of the "granular 
temperature" T in the cross-section yz plane. Height 
profiles of T, averaged on the whole width for each z, 
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Figure 15. (Color online). Vertical profile of the velocity 
fluctuations, averaged over y. 
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Figure 16. (Color online). Transverse profile of the velocity 
fluctuations at the flat base, smoothed over ±2.51) in y. 



are shown in Fig. [15] computed in the bulk of the flow 
above zq — 1.5D. The strong velocity gradient at the 
transition between the bulk and the basal rolling layer 
(see Fig. 13 1 prevents a meaningful computation of T in 
that transition. Fig. [l5]also shows as separate points the 
basal "temperature" value that takes into account only 
contributions from the rolling layer. The main bulk of 
the flow thus rests on top of an effective base with higher 
"granular temperature" and large gradient, which then 
decreases according to a height profile compatible with 
these found in the core flow on bumpy bases, in numerical 
simulations [T5] of thick flows. In these simulations |12| . 
the influence of a bumpy base extends to z — 5D at 
which point the "granular temperature" is maximal (Fig. 
6 of [H]), and then it decreases with the height. Note 
that the temperature profiles reported in Fig. 15 are also 
compatible with those predicted by the kinetic theory |131 
150] and with the early 2D simulation results in |25| . 

The granular temperature T is nearly constant over 
the bulk (Fig. 15 ) in the unidirectional regime 
Fig 



In the 

17] shows a z profile inversion 



convective regime 
between the values at walls and the center that resembles 
Figs. 15a and 15b of [33], where similar profiles were 
computed on a bumpy base with flat frictional walls. The 
highest temperatures occur near the center of the base 
(Fig. 16), the temperature gradient is large near the 
side walls (at least at the base). All these features are 
thus coherent with the idea of a basal layer producing an 
effective bumpy base. 



D. Analogy with Rayleigh-Benard convection 

Matching experimental evidence for secondary flows 
was first seen in |llj . on a bumpy base. Spontaneous 
generation of longitudinal vortices in rapid granular flows 



down rough inclined planes are also been reported in 
|19j . The dense and faster troughs correspond to the 
downward part of the flow, while the dilute and slower 
crest correspond to the upward part. In order to ex- 
plain them an analogy with Rayleigh-Benard convection 
was proposed [13] for granular flows based on consider- 
ations from the kinetic theory for granular gases [6]. A 
three-dimensional linear stability analysis of SFD flows 
reveals that in a wide range of parameters, they are un- 
stable under transverse perturbations. The structure of 
the unstable modes is globally in good agreement with 
the rolls we observe in the main plug of our flows on a 
flat base, despite packing fractions reaching high values 
Vmax > 0.4 in the core in our case (see Fig. 6]), unlike 
the experiments on a rough base reported in 13 where 
Vmax ~ 0.2. These values of v max are similar to those 
obtained in the numerical simulations of [TU] • In [TU] two 
different regimes of stripes are described, but the gran- 
ular temperature is not available. The "dilute" regime 
corresponds to the regime described in [T31 [H?] where the 
dense fast region with downwards motion corresponds to 
a height minimum, while in the "dense" regime it cor- 
responds to a height maximum. The dense regime is 
observed for an average packing fraction V comprised be- 
tween 0.36 and 0.57, while 0.12 < v < 0.42 in the dilute 
regime. It is difficult in our case to know which regime 
correspond our rolls correspond to. The average density 
is in the common range, and the curvature of the surface 
is not a clear indication as the walls could deform it. 

Fig. [17] shows the color-coded map of T and the veloc- 
ity field in the cross-section yz plane. We can see that the 
motion in the bulk part of the transverse plane consists 
of a pair of counter-rotating vortices - The Fig. 5b in flT] 
shows a similar roll orientation. The material moving to- 
wards the base, in the central part is flowing faster in x 
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Figure 17. (Color online). Vector field of the velocities in the transverse yz plane on top of the color/gray-coded "granular 
temperature" T. The data, obtained in the steady and fully developed regime, have been averaged over time and over the 
periodic cell in the direction of the flow. 



direction than the grains rising up on the sides. The av- 
erage density is higher where the flow is downwards and 
smaller where the flow is upwards. We also observe the 
temperature vertical profile inversion reported in Fig. 8d 
of [13]: the temperature gradient is opposed to the trans- 
verse velocity in the downward and upward parts of the 
vortices. 

Interestingly, the Rayleigh-Benard regime is similar to 
the convection that occurs when a granular bed on a 
bumpy base is shaken at high intensity |45) . In such a 
system the shaking and the bumpiness of the base lead to 
a higher granular temperature in the vicinity of the base. 
The granular bed is then heated from below and cooled 



from above. In our system, as shown above, the granu- 
lar layer in contact with the flat frictional base can be 
considered as a bumpy sliding base atop which a sheared 
flow occurs. 



VI. VISCOPLASTIC RHEOLOGY 

A viscoplastic rheology for incompressible flows was 
proposed in [7] , as a 3D extension of the proposal in [3] . 
We expect it to hold in the unidirectional flows case, for 
which |3] was proposed. In j!0| . Borzsonyi et al. have 
shown that the viscoplastic rheology does not hold lo- 
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Figure 18. (Color online). Vertical profile of the inertial num- 
ber I. 



cally for granular convection in the bumpy boundaries 
case. They then propose an extension for this rheology 
to compressible flows. 

This section investigates the situation for the flat fric- 
tional scenario, with an effective basal layer, together 
with a form of the rheology using averaged quantities for 
a global analysis. We average the constitutive equation 
proposed in (7 to match the vertical profile of P(z) pre- 
sented in the previous section, leading to the definition of 
an inertial number I(z) = \{ A f) (z)\ D / \J 'P(z) / p g . In this 
expression the strain rate tensor 7 is averaged at each z 
location over L and W, (7) (z) = 
and the norm used is the same as in [7_ 



>v= 
(i.e. 



\T 



a 2 - 



which recovers the ID expression of I from 



[3] when the 7 tensor is strongly dominated by )■ We 
use the definition of I from [7] in which the grain density 
p g is used for normalization, while the original proposal 
[J used the density of the continuum p c = p g v. Liter- 
ature on this topic shows that both approaches are in 
use (e.g. (24] [46] use p g , [10] use p c ). The averaging 
we propose cancels the transverse v y and vertical v z ve- 
locity components thanks to the symmetry of the inner 
rolls (see Fig. 17 1, recovering an expression that effec- 



tively behaves as for a unidirectional flow without inter- 
nal structure. Moreover, we checked that the norm of 
(7) (z) differs from the norm of its deviator by less than 
0.35%, hence the condition for an incompressible flow is 
satisfied on average (i.e. local dilations that may occur 
along y in the rolls, if any, cancel in any given z slice). 
In these conditions, we expect the constitutive equation 
using the average (I) to hold quite well, which is indeed 
the case (see below). 

Fig. 18 shows the vertical profile of I(z). Above the 



basal layer the average (I(z)) 6 § >z /d>z q /d=i 5 ^ s defined 
in the main bulk of the flow. Oscillations over that av- 
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Figure 19. Fit quality of the theoretical Bagnold profile for 
each angle. 
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Figure 20. (Color online), (a) Effective friction p and (b) 
average packing fraction v vs average inertial number (I) in 
the bulk of the flow. Overlaid thin black lines in (a) are the 
empirical fits mentioned in the main text. 



eraged value are present in the unidirectional regimes, 
matching the layered structure previously mentioned in 
Section [V A[ while for the convective regimes the vertical 
profile does not vary as much I(z) ps (I). Note that for 
the dense flows on bumpy base case, I is assumed to be 
constant on the whole height (Section 8.4.1 of [1]). 

If the constitutive equation of [7] holds in the averaged 
form we propose, we shall recover a velocity profile in the 
form of a Bagnold scaling (eq. 25 of [3]), with A(6) = 
I (I) \J Dcosd matching the constant fitted in the previous 
section. Fig. 



19 shows that this is indeed the case, up to 



a worst-case 5% accuracy. 

From the momentum balance for the main bulk flowing 
on the basal layer, including wall effects, it is possible to 
obtain p((I)), the effective friction coefficient of the main 
bulk on the basal layer: p = tan(9) — pw [Hp — zq) /W 
El 123] • When (I) is plotted against p((I)) the sep- 
aration between the unidirectional and the convective 
regimes is clearly apparent as a discontinuity, see Fig. 
[20] i. The best fit parameters for the constitutive equa- 
tion p((I)) = p s + {p 2 - p s ) I (1 + h/ (I)) proposed in [7J 
(see Fig. [20k) are p s = 0.0046, p 2 = 0.479 and 7 = 0.13 
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in the convective regime, while fi s — 0.16, fi2 = 1-02 and 
Iq = 2.06 in the unidirectional regime, showing that the 
empirical constitutive equation /Lt((7}) changes during the 
transition. The break is similarly visible on the v vs (I) 
profile in Fig. [20] d. Both branches decrease nearly lin- 
early, compatible with Fig. 2 in [46| where the model 
coefficient of restitution and spring stiffness are varied in 
a 2D simulation, and unlike Fig. 4e of [TO] where the 
D(I) relation is built locally and not in averaged form. 

The number I can also be interpreted as the ratio of 
a macroscopic rearrangement time scale over a shearing 
time scale [I]. The observed drastic reduction in I at 
the transition between the unidirectional and convective 
regimes reflects the fact that granular convection rear- 
ranges the grains much faster than slow diffusion within 



20 



the ordered layering. The (I) value at 9 = 18° in Fig 
is consistent with the convective regime despite compu 
tations being performed in the oscillating state, leading 
to the hypothesis that the oscillations are related to the 
onset of convection. That hypothesis will be investigated 
in a future work. 



VII. CONCLUSION 

Our numerical simulations with side walls generate 
SFD flows comparable to the experimental setup [5] with 
a compatible range of angles, distances of establishment 
and velocity profiles. We confirm that the influence of 
the frict ion on the lateral walls is negligible ([2_ and Sec- 
tion 



VI I, but also that walls manifest in other ways a 



long-range influence within the flow ([5] and Fig. 12 1. In 



any case, side walls cannot be ignored even when they are 
far away, especially since channeled flows can be directly 
compared to experiments. Building on these results we 
extrapolate the simulation to larger inclination angles 



and find that distances for reaching the steady states 
exceed the experimental chute length. These regimes 
also correspond to the presence of granular convection, 
whereby grains are circulated within the whole flow, un- 
like the unidirectional regimes where grains mostly re- 
main in a "crystallized" layered structure. 

Compared to the well-studied bumpy base scenario, 
flows on flat frictional surfaces involve a much faster over- 
all velocity, thanks to the presence of a basal layer of 
rolling grains, upon which slides the main bulk of the 
flow. We then interpret that bottommost layer of grains 
as an effective base for the flow bulk and we show that 
in these conditions, the bulk follows a conventional Bag- 
nold scaling. The analogy with an effective rough base 
extends to the presence of a convective regime with sim- 
ilar velocity and density profiles. However, due to the 
increased overall velocity, and owing to the effective base 
being less rigid than a fixed bumpy one, the convection 
rolls appear for lower angles and mass holdups in the flat 
frictional case than in the bumpy one. 

As for the bumpy case, we find that over the effective 
base the bulk of the flow follows on average a viscoplastic 
rheology [7], for each of the SFD regimes. The transition 
between these regimes corresponds to a break in the fric- 
tion n versus inertial number / relation (Fig. 20 1, with 
a drastic reduction in / that matches the effect of the 
secondary rolls (faster grain rearrangement). 

Channeled flows down flat frictional surfaces are well 
adapted for testing granular rheologies numerically and 
studying boundary conditions. 
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